Note: Some of the figures in this document use grid.arrange or ggarrange. The heights/widths arguments in grid.arrange() and ggarrange() may have to be tweaked to get the proportions of the panels to look less lopsided.
Create design matrix using all permunations of sample codes. 1 and 0 corresponds to which side of the comparison the sample in the sample_names column is. The sample names will just correspond to the rownames of the phenotype table. this gets combined in a for-loop later.
This design matrix is referred to as the “pooled comparison strategy” used in the analysis, as it allows us to look at the differences between different regions without resorting entirely to 1v1 comparisons. This way, we can identify gene signatures unique to a particular region compared to the rest, or to identify gene signatures that may be associated with samples containing similar clonal architectures from the previous clonevol analysis.
d1 <- ((combinat::permn(as.character(c('1', '1', '0', '0')))))
d1 <- t( unique( data.frame( matrix(unlist(d1), nrow = factorial(4), byrow = T))))
class(d1) <- "numeric"
d2 <- ((combinat::permn(as.character(c('1', '1', '1', '0')))))
d2 <- t( unique( data.frame( matrix(unlist(d2), nrow = factorial(4), byrow = T))))
class(d2) <- "numeric"
design <- cbind(d1, d2)
colnames(design) <- NULL
rownames(design) <- NULL
# clean up intermediaries
rm(d1, d2)comparison_list <- NULL # store the comparisons here.
deg_list <- NULL # store the results of the comparison here.
# loop through the combinations in the design matrix to do permuted DEG tests in pooled
# comparisons using DESeq2.
for (i in 1:ncol(design)){
sample_code <- design[, i]
ph <- data.frame(sample_name=rownames(case.phenotype), sample_code=sample_code)
dds <- DESeqDataSetFromMatrix(countData = gcounts, colData = ph, design = ~sample_code)
dds <- DESeq(dds)
res <- results(dds)
resOrdered <- as.data.frame(res[order(res$padj),])
resOrdered$ensid <- rownames(resOrdered)
deg_list[[i]] <- resOrdered
comparison_list[[i]] <- ph
}
# clean up intermediaries
rm(resOrdered, dds, res, ph, sample_code)# now take the comparison list into something easier to coerce into a sample v sample string
# for the volcano plot labels
comparison_names <- NULL
for (i in 1:length(comparison_list)){
comparison <- comparison_list[[i]]
s0 <- paste(as.character(comparison[comparison$sample_code == 0,]$sample_name), collapse = ', ')
s1 <- paste(as.character(comparison[comparison$sample_code == 1,]$sample_name), collapse=', ')
deg_list[[i]]$sample_comparison <- as.character(paste(s0, s1, sep=' - vs - '))
comparison_names[i] <- as.character(paste(s0, s1, sep=' - vs - '))
}
names(deg_list) <- comparison_names
rm(s0, s1)
# get top 25 p values and top/bottom logFC values
top_genes <- lapply(deg_list, function(x){
x <- x %>% filter(!is.na(padj)) %>%
filter(padj < 0.5) %>%
left_join(genes, by=c('ensid'))
x1 <- x %>% top_n(15, wt=as.numeric(-log(padj)))
x2 <- x %>% filter(log2FoldChange > 1) %>% top_n(10, wt=abs(log2FoldChange))
x3 <- x %>% filter(log2FoldChange < -1) %>% top_n(10, wt=abs(log2FoldChange))
x <- unique(rbind(x1, x2, x3))
rm(x1, x2, x3)
return(x)
})
top_genes <- do.call(rbind, top_genes)# coerce results to data frame and remove mirrored comparisons
mirrored_comparisons <- c("C, P - vs - A, B", "B, C - vs - A, P", "B, P - vs - A, C")
deg_df <- do.call(rbind, deg_list)
deg_df <- deg_df %>% as_tibble() %>%
filter(!is.na(padj)) %>%
filter(!sample_comparison %in% mirrored_comparisons) %>%
left_join(genes, by=c('ensid'))
top_genes <- top_genes %>%
filter(!is.na(padj)) %>%
filter(!sample_comparison %in% mirrored_comparisons)
summary(as.factor(top_genes$sample_comparison))## A - vs - B, C, P A, B - vs - C, P A, C - vs - B, P A, P - vs - B, C
## 1 33 17 1
## B - vs - A, C, P C - vs - A, B, P P - vs - A, B, C
## 35 32 34
volcano_plot_filt <- ggplot(data=deg_df[! deg_df$sample_comparison %in% c('A - vs - B, C, P', 'A, P - vs - B, C'),]) + aes(x=log2FoldChange, y=-log(padj), color=biotype) +
geom_point() +
geom_label_repel(data=top_genes[! top_genes$sample_comparison %in% c('A - vs - B, C, P', 'A, P - vs - B, C'),], aes(label=gene), show.legend=FALSE, force=3, size=4) +
facet_wrap(~ sample_comparison, scales='free',nrow=1) +
theme_classic() + theme
volcano_plot_filtFirst coerce the deg results from DEseq2 into a list of ranks for each comparison using the method described here: https://stephenturner.github.io/deseq-to-fgsea/
Then, we will use fgsea to perform gene set enrichment on all the MSigDB oncogenic signatures.
deg_ranks <- lapply(deg_list, function(x){
x %>%
left_join(genes, by=c('ensid')) %>%
na.omit() %>%
distinct() %>%
group_by(gene) %>%
summarise(stat=mean(stat)) %>%
deframe()
} )
deg_ranks <- deg_ranks[! names(deg_ranks) %in% mirrored_comparisons]
deg_ranks <- deg_ranks[! names(deg_ranks) %in% c("A - vs - B, C, P","A, P - vs - B, C")]
load('../data/msigdb_gene_signatures.rdata')Generate the enrichment scores & coerce to dataframe
# compute gene set enrichments
all_fgsea <- list(hallmark_fgsea = gsea_lister(rank_list = deg_ranks, sig_list = hallmark_gsigs, sig_name='HALLMARK'),
kegg_fgsea = gsea_lister(rank_list = deg_ranks, sig_list = kegg_gsigs, sig_name = 'KEGG'),
reactome_fgsea = gsea_lister(rank_list = deg_ranks, sig_list = reactome_gsigs, sig_name = 'REACTOME'),
biocarta_fgsea = gsea_lister(rank_list = deg_ranks, sig_list = biocarta_gsigs, sig_name = 'BIOCARTA')
)Select the top pathways, filter irrelevant terms, and order the pathways via hierarchical clustering.
# select top pathways & order pathways via hclust
data_list <- lapply(all_fgsea, function(x){
sig_name <- levels(as.factor(x$signature_name))
x$pathway <- gsub(paste0(sig_name, '_'), '', x=x$pathway)
x$sample_comparison <- factor(x$sample_comparison, levels=names(deg_ranks))
df <- x %>%
filter(padj < 0.05) %>% arrange(padj) %>% group_by(sample_comparison) %>%
top_n(10, wt=-padj) %>% top_n(10, wt=abs(NES)) %>%
arrange(abs(NES), -padj)
tmp <- df %>% dplyr::select(sample_comparison, pathway, NES) %>%
group_by(sample_comparison) %>% tidyr::spread(pathway, NES) %>%
as.data.frame()
names <- tmp$sample_comparison
tmp$sample_comparison <- NULL
tmp <- as.data.frame(lapply(tmp, function(x){
x[is.na(x)] <- 0
return(x) }))
rownames(tmp) <- names
order <- hclust( dist(t(tmp), method = "euclidean"), method = "ward.D" )$order
df$pathway <- factor(df$pathway,levels=colnames(tmp)[order])
df <- df[!is.na(df$pathway),]
df <- df[df$pathway != 'NA',]
rm(names, tmp)
return(df)
})
terms <- c('INFLUENZA', 'CELL_CYCLE', 'MITOTIC', 'S_PHASE', 'G1_S_TRANSITION', 'ASTHMA', 'AUTOIMMUNE_THYROID_DISEASE',
'INFECTION','DISEASE', 'INTESTINAL_IMMUNE_NETWORK_FOR_IGA_PRODUCTION','METABOLISM', 'APICAL_SURFACE','GENESIS',
'43S', 'SYNAP','NEUROTRANSMITTER', 'GABA', 'GLUTAMATE', 'POTASSIUM_CHANNELS', 'CHANNEL_TRANSPORT', 'CELL_LINEAGE',
'LUPUS', 'NEURONAL', 'RP_DEPENDENT_COTRANSLATIONAL_PROTEIN_TARGETING_TO_MEMBRANE', 'TRANSLATION', 'proteasome', 'complement', 'homologous_recombination',
"peptide_chain",'second_messenger', 'REACTIVE_OXYGEN_SPECIES', 'NONSENSE_MEDIATED_DECAY', 'G2M_checkpoint', 'G2_M_CHECKPOINT', 'IMMUNOREGULATORY_INTERACTIONS_BETWEEN',
'RIBOSOME', 'TIGHT_JUNCTION' )
data_filt <- lapply(data_list, function(x){
y <- x
for (term in terms){
y <- y %>% filter(!grepl(term, x=pathway, ignore.case=TRUE))
}
df <- y
tmp <- df %>% dplyr::select(sample_comparison, pathway, NES) %>%
group_by(sample_comparison) %>% tidyr::spread(pathway, NES) %>%
as.data.frame()
names <- tmp$sample_comparison
tmp$sample_comparison <- NULL
tmp <- as.data.frame(lapply(tmp, function(x){
x[is.na(x)] <- 0
return(x) }))
rownames(tmp) <- names
order <- hclust( dist(t(tmp), method = "euclidean"), method = "ward.D" )$order
df$pathway <- factor(df$pathway,levels=colnames(tmp)[order])
df <- df[!is.na(df$pathway),]
df <- df[df$pathway != 'NA',]
rm(names, tmp)
return(df)
})Make one plot per gene set database and visualize as a set of stacked heatmaps.
# set up themes for ggplotting
theme <- theme(text = element_text(size=16),
axis.text.y = element_text(size=18, face='bold'),
axis.text.x = element_text(size=18, face='bold'),
axis.title = element_text(size=24, face = "bold"),
legend.text= element_text(size=18),
legend.title= element_text(size=24, face='bold'),
strip.text=element_text(size=24, face='bold'),
plot.title=element_text(size=28, face='bold'))plot_list <- lapply(data_filt, function(x){
plot <- ggplot(data=x) +
aes(x=sample_comparison, y=pathway, fill=as.numeric(NES)) +
geom_tile(color='white', size=1) +
theme_classic() +
theme + theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank()) +
scale_x_discrete(drop=FALSE) +
labs(x ='', y=levels(as.factor(x$signature_name))) +
scale_fill_gradient2(low = "lightseagreen", mid='white', high = "tomato", name="NES") +
guides(fill = guide_colourbar(barwidth = 1.5, barheight = 25))
return(plot)
})
plot_list[[4]] <- plot_list[[4]] + theme_classic() + theme + rotate_x_text(angle=45)
names(plot_list) <- c('HALLMARK', 'KEGG', 'REACTOME', 'BIOCARTA')
pathway_hm <- ggarrange(plotlist=plot_list, nrow=length(all_fgsea), align='v', common.legend=TRUE, legend='right', heights=c(1.2,0.6,0.5,1.6))surv_pca <- ggplot(data = case_int_pcs) +
aes(x=Dim.1, y=Dim.2, color=log(MonthsOverallSurvival), label=label, shape=VitalStatusFactor) +
geom_density2d(color='black',alpha=0.5, bins=50, aes(shape=NULL)) +
geom_point(alpha=0.8, size=6) +
geom_label_repel(size=8, show.legend=FALSE, force=10, point.padding = unit(0.5, 'lines'), fontface='bold') +
labs(x='PC-1: 14.3% of variance explained', y='PC-2: 9% of variance explained') +
theme_classic() +
theme(plot.margin=margin(0.5,0.5,0.5,0.5, "in")) +
theme + scale_color_gradient(low='ivory', high='black', limits=c(-2,5), breaks=seq(-2,5,length.out = 5)) +
guides(color = guide_colourbar(barwidth = 1.5, barheight = 15)) ## Warning: Removed 3 rows containing missing values (geom_point).
## Warning: Removed 151 rows containing missing values (geom_label_repel).
Objective: get a measure (a p-value) of whether the variance between my four samples is greater than if you were to take the variance four random TCGA samples. To do this, do resamplings of 10% of the total possible combinations of 4 tcga-GBM samples from a pool of 155 samples.
checkpoint.genes <- c("PD1","PDL1","PDL2","TIM3","CD28","CD3","CD4","FOXP3","DCLAMP")
val_cols <- checkpoint.genes
chk <- ind
chk$submitter_id <- as.character(chk$submitter_id)# separate out samples from case study
ref_df <- chk[1:4,]
ref_df <- data.frame(var=sapply(ref_df[, val_cols], var),
samp_type=as.character('test'),
val=val_cols)
set.seed(0)
length(combn(chk[- c(1:4),]$submitter_id, 4, simplify=FALSE))
boot_pool <- chk[- c(1:4),]
# grab all unique combinations of 4 from the pool of 155 TCGA samples
# and then only take 10% of that cause there's too many
combinations <- combn(boot_pool$submitter_id, 4, simplify=FALSE)
combinations <- base::sample(combinations, length(combinations)*0.1)# function to take variance of each combination of 4 for each gene/column
boot_fx <- function(x){
# grab the combination of 4
x <- boot_pool %>%
filter(submitter_id %in% x) %>%
dplyr::select(val_cols)
# take the variance for each gene
x <- data.frame(var=sapply(x, var),
val=val_cols,
samp_type=as.character('random_boot'))
return(x)
}
# do the resampling; using 3 cores (but the more the merrier!), and then bind the rows.
boot_pool <- bind_rows(pbmclapply(combinations, boot_fx, mc.cores=2))
all_boots <- bind_rows(ref_df, boot_pool)Get a measure of whether the variance between my four samples is greater than if you were to take four random TCGA samples.
pvals <- lapply(levels(factor(all_boots$val)), function(i){
x <- all_boots %>% filter(val == i) %>% as.data.frame()
test_point <- x %>% filter(samp_type == "test") %>% dplyr::select(var) %>% deframe()
z <- stats::density(x$var, n=length(x$var))
xval <- as.numeric(as.matrix(as.numeric(z$x)))
yval <- as.numeric(as.matrix(as.numeric(z$y)))
# test right-sided
rpval <- auc(x=xval, y =yval, from=test_point, to=max(xval) )
# test left-sided
lpval <- auc(x=xval, y =yval, from=min(xval), to=test_point )
# pick the smallest of the two left-sided and right-sided tests
# to get the auc (p-value) for the closest end of the distribution
result <- min(c(rpval, lpval))
return(result)
})
names(pvals) <- levels(factor(all_boots$val))
pvals <- unlist(pvals)
pvaldf <- data.frame(pval=pvals,
gene=names(pvals))ind2 <- tidyr::gather(chk, gene, cpm, val_cols)
ind2$gene <- as.factor(ind2$gene)
ind2$cpm <- as.numeric(ind2$cpm)
ind2 <- ind2 %>% inner_join(pvaldf, by=c('gene'))
ind2$sample <- as.factor(ind2$sample)Note: we saved a pre-computed version of this to reduce runtime and will use it in this RMarkdown document.
checkpoint_violins <- ggplot(data=ind2) + aes(y=cpm, x=gene, color=sample) +
#geom_violin(alpha=0.05, aes(color=NULL), fill="#7f7f7f", show.legend = FALSE) +
#geom_boxplot(width=0.25, aes(color=NULL, fill=NULL), show.legend = FALSE ) +
geom_beeswarm(data=ind2[ind2$sample =='TCGA-GBM',], alpha=0.4, size=2, show.legend=FALSE) +
geom_point(data=ind2[ind2$sample !='TCGA-GBM',], size=6, alpha=0.75) +
theme_classic() + theme + labs(y='Normalized Expression', x='', color='Sample') +
# don't add duplicate labels
geom_text(data=ind2[(ind2$pval<0.05) & (ind2$sample=='Primary'),],
y=-5,size=6, show.legend=FALSE,
aes(x=gene, label=paste('P =',round(pval, digits=4)), color=NULL)) +
scale_color_manual(values = cbpal, aesthetics = 'color')Supplementary Figure 3 - overlay onto the PCA
gene_pca <- ggplot(data = ind2) +
aes(x=Dim.1, y=Dim.2, color=cpm, label=label) +
geom_density2d(color="#7f7f7f", bins=50, alpha=0.7) +
geom_point(alpha=0.7, size=5) +
geom_label_repel(size=8, show.legend=FALSE, force=10, min.segment.length = 10, segment.size = 1) +
labs(x='PC-1: 14.3% of variance explained', y='PC-2: 9% of variance explained') +
facet_wrap(~ gene) +
theme_classic() +
theme + scale_color_gradient(low="ivory", high="black", name='Norm. Exprs.', limits=c(-5,10), breaks=seq(-5, 10, by=5)) +
guides(color = guide_colourbar(barwidth = 1.5, barheight = 25)) ## Warning: Removed 3 rows containing missing values (geom_point).
## Warning: Removed 151 rows containing missing values (geom_label_repel).
#pdf('../../figures/fig4.pdf', height=18, width=34)
grid.arrange(pathway_hm, tcga_integration2, ncol=2, widths=c(0.7, 1.5))surv_dens <- ggplot(data = case_int_pcs) +
aes(x=MonthsOverallSurvival) +
geom_density(fill='gray', alpha=0.7, show.legend = FALSE) +
theme_classic() +
theme +
geom_vline(data=case_int_pcs[!is.na(case_int_pcs$label),], size=1.5, color='blue',
aes(xintercept=MonthsOverallSurvival), show.legend = FALSE, linetype='dashed') +
geom_text(data=case_int_pcs[1,],
aes(x = MonthsOverallSurvival, y=0.02), label = 'Patient Survival', color='blue',
show.legend = FALSE,angle=90, size=8, nudge_x =3, fontface='bold') +
scale_y_continuous(limits = c(0,0.05)) ## Warning: Removed 4 rows containing non-finite values (stat_density).
The patient survived almost 1 and a half times longer than the third quartile of tcga survival.
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.16 5.08 11.84 14.15 17.75 88.14 4
Use ECDF on the TCGA cohort to find what the probability is that our patient survived this long by chance.
## [1] "The probability that the patient would have survived 25 months by chance alone is 15.89%"
Taking the pre-treatment and post-treatment responder and non-responder signatures in figure 3b, and ext data fig 7, do single sample GSEA on the relevant gene sets to see if there is any association.
set.seed(0)
zhao_enrichments <- list(pretx_responder_down = gsea_lister(rank_list = deg_ranks, sig_list =zhao_responder_pretx_dn,
sig_name='PRE-TX_DN'),
pretx_responder_up= gsea_lister(rank_list = deg_ranks, sig_list = zhao_responder_pretx_up, sig_name = 'PRE-TX_UP'),
posttx_responder_down = gsea_lister(rank_list = deg_ranks, sig_list = zhao_responder_posttx_dn, sig_name = 'POST-TX_DN'),
posttx_responder_up = gsea_lister(rank_list = deg_ranks, sig_list = zhao_responder_posttx_dn, sig_name = 'POST-TX_UP')
)## # A tibble: 26 x 7
## pathway pval padj NES size sample_comparis… signature_name
## <chr> <dbl> <dbl> <dbl> <int> <chr> <chr>
## 1 XU_GH1_EXOGENOUS_… 0.00149 0.0104 1.80 92 A, B - vs - C, P PRE-TX_DN
## 2 GSE21678_WT_VS_FO… 0.00140 0.0104 1.54 162 A, B - vs - C, P PRE-TX_DN
## 3 GSE17721_LPS_VS_P… 0.00121 0.0169 1.62 180 P - vs - A, B, C PRE-TX_DN
## 4 FAELT_B_CLL_WITH_… 0.0168 0.0454 -1.57 48 C - vs - A, B, P PRE-TX_DN
## 5 MARCHINI_TRABECTE… 0.0235 0.0454 -1.67 18 C - vs - A, B, P PRE-TX_DN
## 6 GSE17721_LPS_VS_P… 0.0264 0.0454 -1.37 184 C - vs - A, B, P PRE-TX_DN
## 7 XU_GH1_EXOGENOUS_… 0.0163 0.0454 -1.52 91 C - vs - A, B, P PRE-TX_DN
## 8 MTOR_UP.N4.V1_DN 0.0213 0.0454 -1.37 155 C - vs - A, B, P PRE-TX_DN
## 9 GSE45739_NRAS_KO_… 0.00256 0.0359 -1.48 174 C - vs - A, B, P PRE-TX_DN
## 10 GSE21678_WT_VS_FO… 0.0292 0.0454 -1.34 160 C - vs - A, B, P PRE-TX_DN
## # … with 16 more rows
znmed_list <- lapply(zhao_enrichments, function(x){
x$sample_comparison <- factor(x$sample_comparison, levels=names(deg_ranks))
df <- x %>% filter(padj < 0.05)
plot <- ggplot(data=df) +
aes(x=sample_comparison, y=pathway, fill=as.numeric(NES)) +
geom_tile(color='white', size=1) +
theme_classic() +
theme + theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank()) +
scale_x_discrete(drop=FALSE) +
labs(x ='', y=levels(as.factor(x$signature_name))) +
scale_fill_gradient2(low = "lightseagreen", mid='white', high = "tomato", name="NES") +
guides(fill = guide_colourbar(barwidth = 1.5, barheight = 25))
return(plot)
})
znmed_list[[4]] <- znmed_list[[4]] + theme_classic() + theme + rotate_x_text(angle=45)
zhao_hm <- ggarrange(plotlist=znmed_list, nrow=length(znmed_list), align='v', common.legend=TRUE, legend='right', heights=c(4,1,1,1.7))#pdf('../../figures/Sfig3_tcga_zhao_ext.pdf', height=18, width=34)
grid.arrange(zhao_hm, gene_pca, ncol=2, widths=c(1, 1.5))## Warning: Removed 1359 rows containing missing values (geom_label_repel).
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Catalina 10.15.1
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] grid parallel stats4 stats graphics grDevices utils
## [8] datasets methods base
##
## other attached packages:
## [1] lattice_0.20-38 ggbeeswarm_0.7.0
## [3] fgsea_1.10.1 Rcpp_1.0.3
## [5] ggpubr_0.2.4 magrittr_1.5
## [7] org.Hs.eg.db_3.8.2 AnnotationDbi_1.46.1
## [9] gridExtra_2.3 ggrepel_0.8.1
## [11] combinat_0.0-8 forcats_0.4.0
## [13] stringr_1.4.0 dplyr_0.8.3
## [15] purrr_0.3.3 readr_1.3.1
## [17] tidyr_1.0.0 tibble_2.1.3
## [19] ggplot2_3.2.1 tidyverse_1.2.1
## [21] edgeR_3.26.8 limma_3.40.6
## [23] DESeq2_1.24.0 SummarizedExperiment_1.14.1
## [25] DelayedArray_0.10.0 BiocParallel_1.18.1
## [27] matrixStats_0.55.0 Biobase_2.44.0
## [29] GenomicRanges_1.36.1 GenomeInfoDb_1.20.0
## [31] IRanges_2.18.3 S4Vectors_0.22.1
## [33] BiocGenerics_0.30.0
##
## loaded via a namespace (and not attached):
## [1] colorspace_1.4-1 ggsignif_0.6.0 htmlTable_1.13.2
## [4] XVector_0.24.0 base64enc_0.1-3 rstudioapi_0.10
## [7] bit64_0.9-7 fansi_0.4.0 lubridate_1.7.4
## [10] xml2_1.2.2 splines_3.6.1 geneplotter_1.62.0
## [13] knitr_1.26 zeallot_0.1.0 Formula_1.2-3
## [16] jsonlite_1.6 broom_0.5.2 annotate_1.62.0
## [19] cluster_2.1.0 compiler_3.6.1 httr_1.4.1
## [22] backports_1.1.5 assertthat_0.2.1 Matrix_1.2-17
## [25] lazyeval_0.2.2 cli_1.1.0 acepack_1.4.1
## [28] htmltools_0.4.0 tools_3.6.1 gtable_0.3.0
## [31] glue_1.3.1 GenomeInfoDbData_1.2.1 fastmatch_1.1-0
## [34] cellranger_1.1.0 vctrs_0.2.0 nlme_3.1-142
## [37] xfun_0.11 rvest_0.3.5 lifecycle_0.1.0
## [40] XML_3.98-1.20 MASS_7.3-51.4 zlibbioc_1.30.0
## [43] scales_1.0.0 hms_0.5.2 RColorBrewer_1.1-2
## [46] yaml_2.2.0 memoise_1.1.0 rpart_4.1-15
## [49] latticeExtra_0.6-28 stringi_1.4.3 RSQLite_2.1.2
## [52] genefilter_1.66.0 checkmate_1.9.4 rlang_0.4.1
## [55] pkgconfig_2.0.3 bitops_1.0-6 evaluate_0.14
## [58] labeling_0.3 htmlwidgets_1.5.1 cowplot_1.0.0
## [61] bit_1.1-14 tidyselect_0.2.5 R6_2.4.1
## [64] generics_0.0.2 Hmisc_4.3-0 DBI_1.0.0
## [67] pillar_1.4.2 haven_2.2.0 foreign_0.8-72
## [70] withr_2.1.2 survival_3.1-7 RCurl_1.95-4.12
## [73] nnet_7.3-12 modelr_0.1.5 crayon_1.3.4
## [76] utf8_1.1.4 rmarkdown_1.17 locfit_1.5-9.1
## [79] readxl_1.3.1 data.table_1.12.6 blob_1.2.0
## [82] digest_0.6.22 xtable_1.8-4 munsell_0.5.0
## [85] beeswarm_0.2.3 vipor_0.4.5